"Programming in Math" updated

The paper on Programming in Mathematical Notation has been updated. (see previous blog entry about Programming in Math)

I've added a page about integration by random sampling. This example adds Greek characters, and a foreign function interface for the random number generator. Also, the translation of several separate sums (defining average, variance, etc) into C++ requires the sums to be merged into a single loop.

[Edit 2/24/2007 - update link]

Cluster algorithm

In the paper, Cluster algorithm for pairwise-interacting particles, Stephen Whitelam and Phillip Geissler describe an algorithm for creating a cluster of particles to move as a group.

The cluster starts with a randomly chosen particle and grows by adding particles with some probability (based on the Boltzmann factor). Particles that strongly interact with cluster are more likely to be part of the cluster. Then a normal trial move is made of the cluster as a whole.

Importantly, the potential and temperature chosen for this operation are not the same as the potential and temperature of the system. In fact, the cluster size can be controlled by adjusting the temperature (at infinite temperature each cluster is only a single particle). This cluster temperature is varied during the simulation run to sample moves with different cluster sizes.

See the paper for more details about controlling the range of energies involved in the cluster (and the cluster interface energies).

Programming in Math

In a previous post on Solving science problems, I commented on the process of scientific computation. One of the steps was converting the solvable equations into a computer program.

As part of making that step more transparent, I've been working on method for Programming in Mathematical Notation (works best with MathML capable browser).

The program takes a content MathML file as input and outputs a presentation MathML file for nicely rendered display in the browser, and a C++ file to compile and run.

[Edit 2/24/2007 - update link to example]

Fast variance optimization

VMC variance optimization usually uses the reweighted variance.

What is reweighting? Consider a set of particle positions obtained from sampling a wavefunction at a given set of parameters. The energy (and the variance of the energy) explicitly depends on the wavefunction parameters by the local energy formula, and implicitly through the particle positions determined by sampling the wavefunction. When optimizing, one would like to vary the parameters without going through all the trouble and expense of sampling new positions. Computing the energy (or variance) change due to the local energy piece is easy - just put the new parameters in and recalculate. Reweighting is the process of correcting the implicit dependence of the energy (or variance) on the wavefunction parameters.


N. D. Drummond and R. J. Needs describe a VMC optimization method in their paper, A variance-minimization scheme for optimizing Jastrow factors, that uses only the explicit dependence of the variance on the parameters (ie, they use the un-reweighted variance). They restrict their scheme to linear parameters in the Jastrow factor, and this makes the resulting parameter surface a quartic. Also, this scheme only requires a single sum over electron positions (the requisite factors can be accumulated during a single run), making the optimization step very fast.

Even though this scheme uses a less accurate approximation for the variance, it appears to give very good results. In fact, the resulting energies for a neon atom are better using the unreweighted variance for optimization than from using the more traditional reweighted variance.

Solving science problems

What is the ideal programming language for MC simulations?

What if we could design a language expressly for MC and QMC simulations - what would it look like? As a first step, let's put programming in the context of solving the whole problem.

Steps to solve a computational science problem

1. Define the problem (ask the question)
   
2. Write down the equations modeling the system under investigation
   
3. Using mathematical transformations and approximations, reduce the equations to a solvable form (numerical analysis)
   
4. Convert into a computer program and run it
   
5. Undertake verification (are steps 2,3, and 4 correct) and validation (are step 1 and the approximations in step 2 correct)
   

Developing techniques for managing each of the steps would be a great help. Current programming languages only deal with step 4. Computer algebra systems (Mathematica, Maple, Maxima, etc) can help with steps 2 and 3, but don't really help in making the jump to a large scale program.

Currently the transformations from steps 2 through 4 are done by hand. I know I make many mistakes in this process. It would be a boon to my productivity if the computer could verify that all the minus signs and factors of 2 are correct.

The verification step often involves understanding the effects of various approximations in step 3. If the program is described as a series of transformations on the equations, hopefully changing the approximation would be a simple change to the equations, and the rest of the process would automatically generate a correct program. This would make it easier to test approximations.

This doesn't say much about the ideal programming language, but it does give some idea of the problem to be solved. I'll have more to say about the ideal programming language later.

Programming Language Notes

On the Sun website, there's an interview with Guy Steele. He talks a little bit about Fortress, which tries to make programs looks more the original problem domain by allowing more mathematical notation. I will be interested to see the results.

Guy Steele is also the principal investigator for the Programming Language Research group. On that page there is a paper titled Object-Oriented Units of Measurement. It's a look the issues involved in making units work in a Java-like OO language.

The last note is an ACM Queue article by Gregory Wilson about Extensible Programming. The article points out that there are advantages to storing the program text in a more structured form, including being able to include diagrams or math notation in the source code.

Whack A Mole

In QMC, a central role is given to the local energy, defined as (1/psi)H psi, where psi is the trial wave function. If psi is an eigenstate, then H psi = E psi by definition, and the local energy becomes a constant. Just as lower energy trial functions are deemed better than higher energy ones, those with fewer fluctutations are better than those with more fluctutations.

One could imagine looking at the fluctuations in the local energy, and using those to determine where best to adjust the trial wave function to reduce the fluctutations (side note: one would want to use psi^2 times the local energy, since fluctuations are weighted by psi^2 in the integrals). This seems similar to the game of Whack-A-Mole, where one attempts to hit mechanical moles that pop up randomly for a brief period of time. Maybe the method could be called Whack-A-Local-Energy (WHALE).

Unfortunately for my naming scheme (and fortunately for everyone else), the method described already exists, and it's called the Energy Fluctutation Potential (EFP) method. It was first described by Stephen Fahy in "Quantum Monte Carlo Methods in Physics and Chemistry" [NATO ASI Ser. C 525 101, 1999; M. P. Nightingale and C. J. Umrigar, eds]. Further papers describing the method:

• Optimal orbitals from energy fluctuations in correlated wave functions, by Fahy and Filippi
  
• Optimization of inhomogeneous electron correlation factors in periodic solids by Prendergast, Bevan, and Fahy
  
• Optimized Jastrow-Slater wave functions for ground and excited states: Application to the lowest states of ethene by Schautz and Filippi
  

The EFP method involves an iteration that starts with a one-body solver (HF or LDA). The resulting orbitals are used in a VMC calculation, and the fluctuations in the local energy form a extra potential that is fed back into the one-body solver to get new orbitals, and the iteration repeats until a self-consistent solution is reached.

(Name note: this paper by Umrigar and Filippi calls it the Effective Fluctuation Potential method. This paper also gives another reference I didn't put on the above list)